#include <stdlib.h>
#include <string.h>
#include <stdio.h>
#include <mpi.h>
#include <math.h>
#include <sys/types.h>
#include <unistd.h>

#define iElemSz 3 
#define ndof 1
#define MaxLineSize 500
#define iFileLineItem 5 
#define dErr pow(10 ,-8) 
#define dLocalErr pow(10 ,-4) 

double **arNodeInfo;
int **arElemInfo;
int **arDomainInfo;
int **arDomainElem;

int iP;
int iConverge;

typedef struct
        {
        double x;
        double y;
        double value;
	int BoundType;
        }stcNode;

typedef struct
        {
	int *arSend;        
	int *arRecv;
        }stcMsg;

stcNode *arPartNode;
int *arPartBound;

stcMsg *arPartMsg;
int *arPartSendNum;
int *arPartRecvNum;
int *arPartElemNum;

int **arPartElem;
int *arIndex;

double **arA;
double *arB;
double **arTA;
double *arTB;
double *arX;

double *arD;		
double *arG;		
double *arTemp;	

int iIterCG, iIterCGSum, iTotalCG;
int iNodeNum;
int iDefCGiter;
char hostname[MaxLineSize];

double **arPartSendValue, **arPartRecvValue;


void AllocPartNode(const int iRank)
	{
	int iNodeNum, iElemNum;

	int i, j, iPartNodeNum=0, iPartBoundNum=0;

	iNodeNum = arNodeInfo[0][0];
	iElemNum = arElemInfo[0][0];
        arPartSendNum = (int *)calloc(iP, sizeof(int));
        arPartRecvNum = (int *)calloc(iP, sizeof(int));

        arIndex = (int *)calloc(iNodeNum+1, sizeof(int));
	arIndex[0] = iNodeNum;

    	for(i=1; i <= iNodeNum ; i++)
		if(arDomainInfo[iRank][i]==1)
			{
			iPartNodeNum++;
			if(arNodeInfo[2][i]==1)
				iPartBoundNum++;
			else
				{
				if(arNodeInfo[3][i]==iRank)
					{
					for(j=0; j<iP; j++)	
						if ((arDomainInfo[j][i]==1)&&(j!=iRank))
							arPartSendNum[j]++;
					}
				else
					arPartRecvNum[(int)arNodeInfo[3][i]]++;
				}
		}
        arPartNode = (stcNode *)calloc(iPartNodeNum+1, sizeof(stcNode));
        arPartBound = (int *)calloc(iPartBoundNum+1, sizeof(int));
	arPartMsg = (stcMsg *)calloc(iP, sizeof(stcMsg));

	arD = (double *)calloc((iPartNodeNum+1),sizeof(double));
	arG = (double *)calloc((iPartNodeNum+1),sizeof(double));
	arTemp = (double *)calloc((iPartNodeNum+1),sizeof(double));

	for(i=0; i<iP; i++)
		{
		arPartMsg[i].arSend = (int *)calloc(arPartSendNum[i]+1, sizeof(int));
		arPartMsg[i].arRecv = (int *)calloc(arPartRecvNum[i]+1, sizeof(int));
		arPartMsg[i].arSend[0]= arPartSendNum[i];
		arPartMsg[i].arRecv[0]= arPartRecvNum[i];

		arPartSendNum[i]=0;
		arPartRecvNum[i]=0;
		}
        arPartNode[0].x = iPartNodeNum;
        arPartBound[0] = iPartBoundNum;

	}


void FillPartNode(const int iRank)
	{
	int iNodeNum, iElemNum;

	int i, j, iPartNodeNum=0, iPartBoundNum=0;

	iNodeNum = arNodeInfo[0][0];
	iElemNum = arElemInfo[0][0];
	
    	for(i=1; i<=iNodeNum ; i++)

		if(arDomainInfo[iRank][i]==1)
			{

			iPartNodeNum++;

			arPartNode[iPartNodeNum].x = arNodeInfo[0][i];	
			arPartNode[iPartNodeNum].y = arNodeInfo[1][i];	
			arIndex[i] = iPartNodeNum;
			if (arNodeInfo[2][i]==1)
				{
				iPartBoundNum++;
				arPartBound[iPartBoundNum] = iPartNodeNum;
				arPartNode[iPartNodeNum].value = arPartNode[iPartNodeNum].x + arPartNode[iPartNodeNum].y ;
				arPartNode[iPartNodeNum].BoundType = 1;	
				}

			else
				{

				if(arNodeInfo[3][i]==iRank)
					{
					for(j=0; j<iP; j++)
						{
						if ((arDomainInfo[j][i]==1)&&(j!=iRank))
							{
							arPartSendNum[j]++;
							arPartMsg[j].arSend[(int)arPartSendNum[j]]=iPartNodeNum;
							}
						}
					}
				else
					{
					arPartRecvNum[(int)arNodeInfo[3][i]]++;
					arPartMsg[(int)arNodeInfo[3][i]].arRecv[(int)arPartRecvNum[(int)arNodeInfo[3][i]]]=iPartNodeNum;
					arPartNode[iPartNodeNum].BoundType = 2;
					}
				}


			}

	}

double *SplitLine(char sLine[], double *arSplitLine)
        {
        #define DLMT " "

        int i = 0;
        char *tokenPtr;
                
        tokenPtr = strtok(sLine, DLMT);
        while (tokenPtr != NULL)
                {
                arSplitLine[i] = atof(tokenPtr);
                tokenPtr = strtok(NULL, DLMT);
                i++;
                }
        arSplitLine[i] = -i;
        return arSplitLine;
        }
           
	
int GetFileInfo(const int iRank, int argc, char **argv)
        {
        int i, j, k;
        int iElemNum, iPNTemp, iPMax = 0;

        FILE *strmNodeInfo, *strmPartInfo, *strmElemInfo;
        char sLine[MaxLineSize];
        char *sNodes;
        char *sParts;

        double arTemp[iFileLineItem] = {0};
        double *arTempPtr = arTemp;

        /* Parse program args to get input file names, and open them */
        if (argc < 3) {
                printf("Usage: %s <nodecount> <partcount>\n", argv[0]);
                return -1;
        }
        sNodes = argv[1];
        sParts = argv[2];

        strcpy(sLine, "./Input/Node");
        strcat(sLine, sNodes);
        if (!(strmNodeInfo = fopen(sLine, "r"))) {
                printf("Cannot open %s for reading.\n", sLine);
                return -1; }

        strcpy(sLine, "./Input/Element");
        strcat(sLine, sNodes);
        if (!(strmElemInfo = fopen(sLine, "r"))) {
                printf("Cannot open %s for reading.\n", sLine);
                return -1; }

        strcpy(sLine, "./Input/Part");
        strcat(sLine, sNodes);
        strcat(sLine, "-");
        strcat(sLine, sParts);
        if (!(strmPartInfo = fopen(sLine, "r"))) {
                printf("Cannot open %s for reading.\n", sLine);
                return -1; }

	fseek(strmNodeInfo, 0, SEEK_SET);
	fseek(strmElemInfo, 0, SEEK_SET);
	fseek(strmPartInfo, 0, SEEK_SET);

	fgets(sLine, MaxLineSize, strmNodeInfo);
        sscanf(sLine,"%d",&iNodeNum);

        arNodeInfo = (double **)calloc(4, sizeof(double*));
        for(i =0; i < 4 ; i++)
                arNodeInfo[i] = (double *)calloc(iNodeNum+1, sizeof(double));

	arNodeInfo[0][0] = iNodeNum;

	for (i=1; i<=iNodeNum; i++)
		{
		fgets(sLine, MaxLineSize, strmNodeInfo);
		arTempPtr = SplitLine(sLine,arTemp);
		for (j=0; j<3; j++)
			arNodeInfo[j][i] = arTempPtr[j+1];		

		fgets(sLine, MaxLineSize, strmPartInfo);
		arNodeInfo[3][i] = atoi(sLine);		
		if (arNodeInfo[3][i] > iPMax)
			iPMax = arNodeInfo[3][i];		
		}

	if (iP != iPMax+1)
		{
		if (iRank==0)
			{
			printf("\n Error! Partition file has %d partitions. It does not match your %d processors requested.\n", iPMax + 1, iP);
			printf(" Program Terminated. \n");
			}
		exit(EXIT_SUCCESS);
		}

	
	arDomainInfo = (int **)calloc(iP, sizeof(int *));
        for(i =0; i < iP ; i++)
		arDomainInfo[i] = (int *)calloc(iNodeNum+1, sizeof(int));
	arDomainInfo[0][0]=iNodeNum;

	fgets(sLine, MaxLineSize, strmElemInfo);
        sscanf(sLine,"%d",&iElemNum);

	arElemInfo = (int **)calloc(iElemSz+1, sizeof(int *));
    	for(i =0; i < iElemSz ; i++)
		arElemInfo[i] = (int *)calloc(iElemNum+1, sizeof(int));

	arElemInfo[0][0]=iElemNum;

	arDomainElem = (int **)calloc(iP, sizeof(int*));
        for(i =0; i < iP ; i++)
		arDomainElem[i] = (int *)calloc(iElemNum+1, sizeof(int));
	arDomainElem[0][0]=iElemNum;

	for (i=1; i<=iElemNum; i++)
		{
		fgets(sLine, MaxLineSize, strmElemInfo);
		arTempPtr = SplitLine(sLine,arTemp);
		for (j=0; j<iElemSz; j++)
			{
			arElemInfo[j][i] = arTemp[j+1];
			iPNTemp = arNodeInfo[3][arElemInfo[j][i]];
			arDomainElem[iPNTemp][i] = 1;
			for(k=1; k<=iElemSz-1; k++)
				arDomainInfo[iPNTemp][(int)arTemp[((j+k)%iElemSz)+1]]=1;
			}
		}
	
	return 0;
	}


void Feeldof(int *node, int *nodeIndex)
        {
        int i, j, k, iTemp;
        k = 0;
        
        for (i=0; i<iElemSz; i++)
                {
                iTemp = (node[i] - 1) * ndof;
                for (j=1; j <= ndof; j++)
                        {
                        nodeIndex[k] = iTemp + j;
                        k++;
                        }
                }
        }  


void LocalMeshToSystem(const int iRank)
        {
	int i, j, k, ii, jj;
	int  iLoc, iElemNum, iPartNodeNum;
	double dDet, dTemp;
	double x[iElemSz], y[iElemSz];
	int node[iElemSz], nodeIndex[iElemSz];
	double arK[iElemSz][iElemSz];

	iElemNum = arElemInfo[0][0];
	iPartNodeNum = (int)arPartNode[0].x;

        arA = (double **)calloc(iPartNodeNum+1, sizeof(double));
        arTA = (double **)calloc(iPartNodeNum+1, sizeof(double));
        for(i =0; i < iPartNodeNum+1 ; i++)
		{
                arA[i] = (double *)calloc(iPartNodeNum+1, sizeof(double));
                arTA[i] = (double *)calloc(iPartNodeNum+1, sizeof(double));
		}
        arB = (double *)calloc(iPartNodeNum+1, sizeof(double));
        arTB = (double *)calloc(iPartNodeNum+1, sizeof(double));
        arX = (double *)calloc(iPartNodeNum+1, sizeof(double));

        for (i=1; i<=iElemNum; i++)
                {
                if (arDomainElem[iRank][i]==1)
			{
                	for (j = 0; j<iElemSz; j++)
				{
				iLoc = arIndex[(int)arElemInfo[j][i]];
				node[j] = iLoc;
				x[j] = arPartNode[iLoc].x;
				y[j] = arPartNode[iLoc].y;
				}

			Feeldof(node, nodeIndex);


        		dDet = 0.5*(((x[1] * y[2]) - (x[2] * y[1])) - (x[0] * (y[2] - y[1])) + (y[0] * (x[2] - x[1])));
        		dTemp = 1 / (4 * dDet);

        		arK[0][0] = dTemp * ((x[2] - x[1])*(x[2] - x[1]) + (y[1] - y[2])*(y[1] - y[2]));
        		arK[0][1] = dTemp * ((x[2] - x[1])*(x[0] - x[2]) + (y[1] - y[2])*(y[2] - y[0]));
        		arK[0][2] = dTemp * ((x[2] - x[1])*(x[1] - x[0]) + (y[1] - y[2])*(y[0] - y[1]));
        
        		arK[1][0] = arK[0][1];
        		arK[1][1] = dTemp * ((x[0] - x[2])*(x[0] - x[2]) + (y[2] - y[0])*(y[2] - y[0]));
        		arK[1][2] = dTemp * ((x[0] - x[2])*(x[1] - x[0]) + (y[2] - y[0])*(y[0] - y[1]));
                
        		arK[2][0] = arK[0][2];  
        		arK[2][1] = arK[1][2];
        		arK[2][2] = dTemp * ((x[1] - x[0])*(x[1] - x[0]) + (y[0] - y[1])*(y[0] - y[1]));

                        for (j=0; j<iElemSz; j++)
                                {
				ii = nodeIndex[j];
				if (arPartNode[ii].BoundType==0)
					{
                                	for (k=0; k<iElemSz; k++)
						{
                        			jj = nodeIndex[k];
                                        	arA[ii][jj] = arA[ii][jj] + arK[j][k];
						}
					}
				else
					{
					arA[ii][ii] = 1;
					arB[ii] = arPartNode[ii].value;
					}
				}  

                        }
		   
		}

        }

//MODIFIED: 21 May 08 - Ax = b -> PAx = Pb
//P is the preconditioning matrix for preconditioned conjugate gradient
void MultTranspose()
	{
		int i, j, k;
		int iPartNodeNum;
		double curPInverse;

		iPartNodeNum = arPartNode[0].x;

		//Precondition the matrix

		//Multiply arA by P inverse (P is the diagnol matrix with the same
		//diagnol as A and 0's elsewhere)
 		for(i=1;i <= iPartNodeNum;i++)
 		{
        	curPInverse = 1/arA[i][i];
        	for(k=1;k<=iPartNodeNum;k++)
        	{
        		arA[i][k] = arA[i][k] * curPInverse; 
        	}
                
        	//Multiply arB by P inverse
        	arB[i] = arB[i] * curPInverse;
		}

		for (i = 1; i <= iPartNodeNum ; i++)
		{
                for (j=1; j<=iPartNodeNum ; j++)
                {
                	arTB[i] = arTB[i] + arA[j][i] * arB[j];
                    for (k=1; k<=iPartNodeNum; k++)
                    {
                        arTA[i][j] = arTA[i][j] + arA[k][i] * arA[k][j];
					}
					
                }
                

    	} 
	}  

void UpdateTBVec()
	{
	int i, j;
	int iPartNodeNum;

	iPartNodeNum = arPartNode[0].x;

	for (i=1; i<=iPartNodeNum; i++)
		{
		if (arPartNode[i].BoundType==2) 
                        arTB[i] = arPartNode[i].value;	
		if (arPartNode[i].BoundType!=1)
			arX[i] = arPartNode[i].value;
		}
	}

double*MultVector(double*arTempVector, const int iNum, double*arTempAnswer)
	{
        int i, j;

        for (i=1; i<=iNum; i++)
                {
                arTempAnswer[i] = 0.0;
                for (j=1; j<=iNum; j++)
                	arTempAnswer[i] = arTempAnswer[i] + arTA[i][j] * arTempVector[j];
                }

        return arTempAnswer;
	}


double InnerProduct(double*arTempVector1, double*arTempVector2, const int iNum)
	{
        int i;
        double dInner = 0.0;

        for (i=1; i<=iNum; i++)
        	dInner = dInner + arTempVector1[i] * arTempVector2[i];

        return dInner;
	}


void CG()        
	{
        int i, j, k,iCG;
	int iPartNodeNum;
        double Temp1, Temp2, Temp3, Temp4, dS;

	iPartNodeNum = arPartNode[0].x;
	iIterCG=0;
	if (iDefCGiter == 0) 
		iCG = iPartNodeNum;
	else
		iCG = iDefCGiter;	

	/*printf (" iCG= %d \n",iCG);*/

        for(i=1; i<=iPartNodeNum; i++)
        	{
                arD[i] = 0.0;
	/*	arX[i] = 0.0;  */
                arG[i] = -arTB[i];
                }

     /*  for(i=1; i<=iPartNodeNum; i++) */  
     /*	  for(i=1; i<=2; i++) */

         for(i=1; i<=iCG; i++)

        	{
                Temp1 = InnerProduct(arG, arG, iPartNodeNum);
		arG = MultVector(arX, iPartNodeNum, arG);

                for(j=1; j<=iPartNodeNum; j++)
                	arG[j] = arG[j] - arTB[j];

                Temp2 = InnerProduct(arG, arG, iPartNodeNum);         
                if (Temp2 < dLocalErr)
                	break;  

                for(j=1; j<=iPartNodeNum; j++)
                	arD[j] = -arG[j] + (Temp2 / Temp1) * arD[j];

                Temp3 = InnerProduct(arD, arG, iPartNodeNum);
		
                Temp4 = InnerProduct(arD,MultVector(arD, iPartNodeNum, arTemp), iPartNodeNum);

                dS = -Temp3 / Temp4;           
		iIterCG++;

                for(j=1; j<=iPartNodeNum; j++)
                	arX[j] = arX[j] + dS * arD[j];
                }
	}

void UpdateValues()
	{
	int i;
	int iPartNodeNum;
	iConverge = 0;
	iPartNodeNum = arPartNode[0].x;
	for(i=1; i<=iPartNodeNum; i++)
		if (arPartNode[i].BoundType==0)
			{	
			iConverge = iConverge + (dErr<fabs(arPartNode[i].value - arX[i]));
			arPartNode[i].value = arX[i];
			}
	}

void CommAloc()
	{
	int i;
	int iSendSz, iRecvSz;

	arPartSendValue = (double **)calloc(iP, sizeof(double)); 
	arPartRecvValue = (double **)calloc(iP, sizeof(double)); 

	for(i=0; i<iP; i++)
		{
		iSendSz = arPartMsg[i].arSend[0];
		iRecvSz = arPartMsg[i].arRecv[0];

		if(iSendSz>0 && iRecvSz>0)
			{
			arPartSendValue[i] = (double *)calloc(iSendSz+1, sizeof(double)); 
			arPartRecvValue[i] = (double *)calloc(iRecvSz+1, sizeof(double)); 
			}
		}
	}


void Communication()
	{
	int i, j;
	int iSendSz, iRecvSz;
	MPI_Status status;

	for(i=0; i<iP; i++)
		{
		iSendSz = arPartMsg[i].arSend[0];
		iRecvSz = arPartMsg[i].arRecv[0];

		if(iSendSz>0 && iRecvSz>0)
			{
			for(j=1; j<=iSendSz; j++)
				arPartSendValue[i][j]=arPartNode[(int)arPartMsg[i].arSend[j]].value;
	
			MPI_Sendrecv(arPartSendValue[i],iSendSz+1, MPI_DOUBLE, i, 123, arPartRecvValue[i],
			iRecvSz+1, MPI_DOUBLE,i , 123, MPI_COMM_WORLD,&status);
		
			for(j=1; j<=iRecvSz; j++)
				arPartNode[(int)arPartMsg[i].arRecv[j]].value = arPartRecvValue[i][j];

			}
		}
	}


void PrintInfo(const int iRank, const int iP, int iStart, int iEnd)
        {
        int i, j;
        int *arCount, *arDisp;
	double *arFinalX, *arPartX, *arFinalSortX;
	int *arPartIndex, *arFinalIndex;
	double dTemp;
	arCount = (int*) calloc(iP, sizeof(int));
	arDisp = (int*) calloc(iP, sizeof(int));
	arFinalX = (double*) calloc(iNodeNum, sizeof(double));
	arFinalSortX = (double*) calloc(iNodeNum, sizeof(double));
	arFinalIndex = (int*) calloc(iNodeNum, sizeof(int));
	
	for(i=1; i<=iNodeNum; i++)
		arCount[(int)arNodeInfo[3][i]]++;


	arPartX = (double *) calloc(arCount[iRank], sizeof(double));
	arPartIndex = (int *) calloc(arCount[iRank], sizeof(int));

	j = 0;
	for(i=1; i<=iNodeNum; i++)
		if(arNodeInfo[3][i]==iRank)
			{
			arPartX[j] = arPartNode[arIndex[i]].value;
			arPartIndex[j] = i;
			j++;
			}

	for(i=0; i<iP-1; i++)
		arDisp[i+1] = arDisp[i] + arCount[i];

        MPI_Gatherv(arPartX, arCount[iRank], MPI_DOUBLE, arFinalX, arCount, arDisp, MPI_DOUBLE, 0, MPI_COMM_WORLD);
        MPI_Gatherv(arPartIndex, arCount[iRank], MPI_INT, arFinalIndex, arCount, arDisp, MPI_INT, 0, MPI_COMM_WORLD);

	for(i=0; i<iNodeNum; i++)
		arFinalSortX[arFinalIndex[i]-1] = arFinalX[i];		

        if (iRank==0)
                {  
                if (iStart<1)
                        iStart = 1;
        
                if (iEnd+1 > iNodeNum)
                        iEnd = iNodeNum;
        
                if (iStart<=iEnd)
                        for(i=iStart-1; i<iEnd; i++)
                        {
                                printf("%f \n", arFinalSortX[i]);
/*                              printf("X[%d] = %f \n", i+1, arFinalX[i]);   */
                        }
                }
        }


void PrintReport(int iP, int iRank, double tStop, double tTotal, double tSumCG, double tCGTotal,int iIterCGSum,int iTotalCG, double 
tSumComm, int iIter)
	{
	int i;
	if(iRank==0)
		{
		printf(" ----------------------------------------------------------------------\n");		
		printf(" Host: %s \n",hostname);
		printf(" #Nodes: %d \n",arIndex[0]);
		printf(" #Elements: %d \n",arElemInfo[0][0]);
		printf(" #Partitions: %d \n",iP);
		printf(" Parallel Time: %f seconds\n", tTotal);
		printf(" CG Time: %f seconds \n", tCGTotal);
		if (iDefCGiter != 0)
                printf(" #CG iterations defined per step: %d \n", iDefCGiter);
		printf(" #CG Total Iterations: %d \n", iTotalCG);
		printf(" #Communications: %d \n", iIter);
		printf(" Tolerance: %E \n", dErr);
		printf(" ----------------------------------------------------------------------\n\n");
		printf(" ----------------------------------------------------------------------\n");
		printf("|  Processor#  |  Total Time(sec)  |  Comm Time(sec)|  Total CG Iter#  |\n");
		printf("|--------------|-------------------|----------------|------------------|\n");
		}

	for(i=0 ;i<iP ;i++)
	        {
		if (i==iRank)
			{
		printf("     P%d             %f             %f         %d  \n",i,tStop,tSumComm,iIterCGSum);
			printf("|--------------|-------------------|----------------|------------------|\n");
			}
		MPI_Barrier(MPI_COMM_WORLD);
		}

	}


void PrintReportOLD(int iP, int iRank, double tStop, double tTotal, double tSumCG, double tCGTotal,int iIterCGSum,int iTotalCG, double
tSumComm, int iIter)
        {
        int i;
/*	FILE *fp1;
	char sLine[MaxLineSize];	
	strcpy(sLine, "./Result/FinalMetis373-4");
        strcat(sLine, iDefCGiter);

  	fp1 = fopen(sLine, "w");
  	fprintf (fp1, "<?\n");
*/
        MPI_Barrier(MPI_COMM_WORLD);

        if(iRank==0)
                {
		printf("\n %d \n", iDefCGiter);
                printf("\n T:%f \n", tTotal);
                printf("%d \n", iIter);
		}

        MPI_Barrier(MPI_COMM_WORLD);

        for(i=0 ;i<iP ;i++)
                {
                if (i==iRank)
                        {
                	printf("%f \n",tSumComm);
			}
                MPI_Barrier(MPI_COMM_WORLD);
                }

	/* fclose(fp1); */

        }



int main(int argc, char **argv)
        {
        int i, iIter=0;
	int iRank, iStop, ret;
	
	double tStartCG, tStopCG, tSumCG=0, tCGTotal;
	double tStartComm, tStopComm, tSumComm=0, tCommTotal;
	double tStart, tStop, tTotal;
	

        MPI_Init(&argc,&argv);
        MPI_Comm_rank(MPI_COMM_WORLD,&iRank);
        MPI_Comm_size(MPI_COMM_WORLD,&iP);


	gethostname(hostname, MaxLineSize); 
	printf("[%d/%d]: %s\n", iRank, iP, hostname);
        if (argc>3)
                {
                if ( strcasecmp(argv[3],"-i") == 0)
                        iDefCGiter = atoi(argv[4]);
		}


        if (argc>6)
                {
                if ( strcasecmp(argv[6],"-i") == 0)
                        iDefCGiter = atoi(argv[7]);
                }

        iIter = 0;

	ret = GetFileInfo(iRank, argc, argv);

        if (ret != 0) 
		{
                printf("Read failed at rank %d.\n", iRank);
                MPI_Abort(MPI_COMM_WORLD, ret);
                return ret;
        	}
		
	AllocPartNode(iRank);

	FillPartNode(iRank);

	LocalMeshToSystem(iRank);

	CommAloc();

	MultTranspose();

	iIterCGSum = 0;

tStart = MPI_Wtime();

	do 
		{
tStartCG = MPI_Wtime();
		CG();
tStopCG = MPI_Wtime() - tStartCG;
		tSumCG = tSumCG + tStopCG;
		iIterCGSum = iIterCGSum + iIterCG;

		UpdateValues();

		MPI_Allreduce( &iConverge, &iStop, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

tStartComm = MPI_Wtime();
		Communication();
tStopComm = MPI_Wtime() - tStartComm;
		tSumComm = tSumComm + tStopComm;
		UpdateTBVec();
		iIter++;
		}while(iStop!=0);

tStop = MPI_Wtime() - tStart;

        if (argc>4)
                {
                if ( strcasecmp(argv[3],"-d") == 0)
                        PrintInfo(iRank, iP, atoi(argv[4]), atoi(argv[5]));
		}

        if (argc>5)
                {
                if ( strcasecmp(argv[5],"-d") == 0)
                       PrintInfo(iRank, iP, atoi(argv[6]), atoi(argv[7]));
                }

	MPI_Reduce( &tStop, &tTotal, 1, MPI_DOUBLE, MPI_MAX,0, MPI_COMM_WORLD);
	MPI_Reduce( &tSumCG, &tCGTotal, 1, MPI_DOUBLE, MPI_MAX,0, MPI_COMM_WORLD);
	MPI_Reduce( &tSumComm, &tCommTotal, 1, MPI_DOUBLE, MPI_MAX,0, MPI_COMM_WORLD);
	MPI_Reduce( &iIterCGSum, &iTotalCG, 1, MPI_INT, MPI_MAX,0, MPI_COMM_WORLD);

	MPI_Barrier(MPI_COMM_WORLD);

	PrintReport(iP, iRank, tStop, tTotal, tSumCG, tCGTotal, iIterCGSum, iTotalCG,tSumComm, iIter);
	//PrintInfo(iRank,iP,1,iNodeNum);
        MPI_Barrier(MPI_COMM_WORLD);

        MPI_Finalize();
        return 0;
        }
